﻿using System;
using System.Collections.Generic;
using System.Data;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
using DevExpress.XtraCharts.Design;
using MathNet.Numerics.LinearAlgebra.Double;

namespace AnalysisSystem
{
    class Smooth
    {
        public static void SmoothMethod(ref DataTable table, string xName, string yName, int point, int times)
        {
            List<double> y = new List<double>();
            Dictionary<int, bool> isNullValue = new Dictionary<int, bool>();
            for(int i=0; i<table.Rows.Count; i++)
            {
                double yValue;
                if (table.Rows[i][yName].Equals(DBNull.Value))
                {
                    isNullValue.Add(i, false);
                    continue;
                }
                else
                {
                    yValue = Convert.ToDouble(table.Rows[i][yName]);
                    if(yValue.Equals(0))
                    {
                        isNullValue.Add(i, false);
                        continue;
                    }
                    else isNullValue.Add(i, true);
                }
                y.Add(yValue);
            }
            List<double> DValue;
            List<double> bValue;
            List<double> yResult = SmoothCalculation(y, point, times, out DValue, out bValue);
            if (table.Columns.Contains("平滑值")) table.Columns.Remove("平滑值");
            table.Columns.Add("平滑值");
            if (table.Columns.Contains("D Value")) table.Columns.Remove("D Value");
            table.Columns.Add("D Value");
            if (table.Columns.Contains("b Value")) table.Columns.Remove("b Value");
            table.Columns.Add("b Value");
            for(int i=0, j=0; i<table.Rows.Count; i++)
            {
                if (isNullValue[i])
                {
                    table.Rows[i]["平滑值"] = yResult[j];
                    table.Rows[i]["D Value"] = DValue[j];
                    table.Rows[i]["b Value"] = bValue[j];
                    j++;
                }
                else
                {
                    table.Rows[i]["平滑值"] = DBNull.Value;
                    table.Rows[i]["D Value"] = DBNull.Value;
                    table.Rows[i]["b Value"] = DBNull.Value;
                }
                //MessageBox.Show(string.Format("原数据：{0}, 平滑值: {1}", table.Rows[i]["日产气（万方/天）"].ToString(), table.Rows[i]["平滑值"]));
            }
        }
        private static List<double> SmoothCalculation(List<double> y, int point, int times, out List<double> D, out List<double> b)
        {
            List<double> smoothY = new List<double>();
            D = new List<double>();
            b = new List<double>();
            for(int i=0; i<point/2; i++)
            {
                smoothY.Add(y[i]);
                //无效值，临时填充使用
                D.Add(65535);
                b.Add(65535);
            }
            for(int i=point/2; i<y.Count-point/2; i++)
            {                
                List<double> dataY = new List<double>();
                for(int j=(-point/2); j<=point/2; j++)
                {                    
                    dataY.Add(y[i + j]);
                }
                double firstDerivative, secondDerivative;
                double value = SGAlgorithm(dataY, point, times, out firstDerivative, out secondDerivative);
                smoothY.Add(value);
                D.Add((-1 / value) * firstDerivative);
                b.Add((secondDerivative * value) / (firstDerivative * firstDerivative) - 1);
            }
            for(int i=y.Count - point/2; i < y.Count; i++)
            {
                smoothY.Add(y[i]);
                //无效值，临时填充使用
                D.Add(65535);
                b.Add(65535);
            }
            return smoothY;
        }

        private static double SGAlgorithm(List<double> ydata, int point, int times, out double firstDerivative, out double secondDerivative)
        {
            List<double> xdata = FillMatrixData(point, times);
            DenseMatrix xmatrix = new DenseMatrix(point, times, xdata.ToArray());
            DenseMatrix ymatrix = new DenseMatrix(point, 1, ydata.ToArray());
            DenseMatrix amatrix = (DenseMatrix)(xmatrix.Transpose().Multiply(xmatrix)).Inverse();
            amatrix = (DenseMatrix)amatrix.Multiply(xmatrix.Transpose()).Multiply(ymatrix);
            firstDerivative = amatrix[1, 0];
            secondDerivative = amatrix[2, 0];
            DenseMatrix yResult = (DenseMatrix)xmatrix.Multiply(amatrix);
            return yResult[point/2, 0];            
        }

        private static List<double> FillMatrixData(int point, int times)
        {
            List<double> matrixData = new List<double>();
            for(int i=0; i<times; i++)
            {
                for(int j=(-point/2); j<=point/2; j++)
                {
                    double value = Math.Pow(j, i);
                    matrixData.Add(value);
                }
            }
            return matrixData;
        }
    }
}
